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Abstract 

Fourier reconstruction algorithms significantly outperform conventional 
back-projection algorithms in terms of computation time. In photoacous- 
tic imaging, these methods require interpolation in the Fourier space do- 
main, which creates artifacts in reconstructed images. We propose a novel 
reconstruction algorithm that applies the one-dimensional nonuniform fast 
Fourier transform to photoacoustic imaging. It is shown theoretically and 
numerically that our algorithm avoids artifacts while preserving the com- 
putational effectiveness of Fourier reconstruction. 

Key— words. Image reconstruction, photoacoustic imaging, planar mea- 
surement geometry, fast algorithm, nonuniform FFT. 

1 Introduction 

Photoacoustic imaging (PAI) is a novel promising tool for visualizing light ab- 
sorbing structures in an optically scattering medium, which carry valuable in- 
formation for medical diagnostics. It is based on the generation of acoustic 
waves by illuminating an object with pulses of non- ionizing electromagnetic ra- 
diation, and combines the high contrast of pure optical and the high resolution 
of ultrasonic imaging. The method has demonstrated great promise for a vari- 
ety of biomedical applications, such as imaging of animals [1, 2], early cancer 
diagnostics [3, 4], and imaging of vasculature [5, 6]. 

When an object is illuminated with short pulses of non- ionizing electromag- 
netic radiation, it absorbs a fraction of energy and heats up. This in turn induces 
acoustic (pressure) waves, that are recorded with acoustic detectors outside of 
the object. Other than in conventional ultrasound imaging, where the source of 
acoustic waves is an external transducer, in PAI the source is the imaged object 
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itself. The frequency bandwidth of the recorded signals is therefore generally 
broad and depends on the size and the shape of illuminated structures. 

1.1 Planar recording geometry 

Throughout this paper we assume a planar recording geometry, where the acous- 
tic signals are recorded with omnidirectional detectors arranged on planes (or 
lines), see Figure 1. The planar geometry is of particular interest since it can 
be realized most easily in practical applications. The recorded acoustic sig- 
nals are then used to reconstruct the initially generated acoustic pressure which 
represents optically absorbing structures of the investigated object. 




Figure 1: Photoacoustic imaging for planar recording geometry. The 

object is illuminated by a pulse of electromagnetic radiation, and reacts with 
an expansion. Induced acoustic waves are measured with an array of acoustic 
detectors arranged on a plane (or a line) and used to from an image of the 
object. 

For the planar recording geometry, two types of theoretically exact recon- 
struction formulas have been reported: Temporal back-projection [7, 8, 9, 10] 
and Fourier domain formulas [11, 7, 9, 12, 13, 14, 15]. Numerical implemen- 
tations of those formulas often lead to fast and accurate image reconstruction 
algorithms. 

In temporal back-projection formulas the signals measured at time t are back 
projected over spheres of radius Vst with the detector position in the centre {vs 
denotes the speed of sound). In Fourier domain formulas this back-projection is 
performed by interpolation in the frequency domain. Reconstruction methods 
based on Fourier domain formulas are attractive since they reconstruct an x 
N X N image in 0{N^ log N) floating point operations by using of the Fast 
Fourier Transform (FFT). Straightforward implementations of back-projection 
type formulas, on the other hand, require 0{N^) operation counts, see [16, 17]. 

The standard FFT algorithm assumes sampling on an equally spaced grid 
and therefore, in order to implement the Fourier domain formulas, interpolation 
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in the Fourier space is required. Interpolation in the Fourier domain is a critical 
issue, and creates artifacts in reconstructed images, see the examples in Section 
5. One obtains significantly better results by increasing the sampling density 
in the Fourier space. This is achieved by either zero-padding [13] or by sym- 
metrizing the recorded signals around t = (which is equivalent to using the 
fast cosine transform instead of the FFT). In this paper, we propose an efficient 
reconstruction algorithm that uses the nonuniform (or unevenly spaced) FFT 
[18, 19, 20, 21, 22, 23] and further increases the quality of reconstruction. 

1.2 Prior work and innovations 

The nonuniform FFT has been applied to a variety of medical imaging prob- 
lems, such as standard X-ray CT, magnetic resonance imaging, and diffraction 
tomography [24, 25, 26], and has also been used implicitly in gridding algo- 
rithms [27, 28]. All those algorithms deal with the problem of recovering a two 
(or higher) dimensional object function from samples of its multidimensional 
Fourier transform on a non-cartesian grid. 

Our approach is conceptually quite different to the above mentioned ref- 
erences: The special structure of our problem allows to perform several one- 
dimensional nonuniform FFTs instead of a single higher dimensional one. This 
leads to a reduced numerical cost, compared to the above algorithms. The pro- 
posed algorithm is more closely related to a reconstruction algorithm for X-ray 
CT suggested in [29, Section 5.2], which also evaluates the Fourier transform on 
irregular samples by means of the one-dimensional FFT. 

1.3 Outline 

This article is organized as follows: In Section 2 we present the mathematical 
basics of Fourier reconstruction in PAL In Section 3 we review the nonuniform 
FFT which is then used to derive the nonuniform FFT based reconstruction 
algorithm in Section 4. In Section 5 we present numerical results of the proposed 
algorithm and compare it with existing Fourier and back projection algorithms. 
The paper concludes with a discussion of some issues related to sampling and 
resolution in the Appendix. 

2 Photoacoustic Imaging 

Let Co^(H) denote the space of smooth functions with bounded support in the 
half space EI := R^~^ x (0, oo), d > 2. Consider the initial value problem 

{d^ - A) p(x, t) = , (x, t) G X (0, oo) , 

p(x,0)=/(x), xGR^ 
dtp{^, 0) = , X G M"^ , 
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with / G Co^(H). Here A denotes the Laplacian in and dt is the derivative 
with respect to t. We write x = (x, x G M^"^, ^ G M, and define the operator 
Q : C^(H) ^ C^(R^) by 



(Q/)(x,t) := 



p(x,^ = 0,t), ift >0, 
, otherwise . 



Photoacoustic imaging for planar recording geometry is concerned with recon- 
structing / G Cq^(EI) from incomplete and possibly erroneous knowledge of Q / . 
Of practical interest are the cases d = 2 and d = 3, see [30, 31, 32, 33]. 

2.1 Exact inversion formula 

The operator Q can be inverted analytically by means of the exact inversion 
formula 

2Ky{¥Qif) (K,,sign{Ky)J\K,\^+KA 
{Ff){K,,Ky) = ^ , ^ ^ (1) 

where {Kx^ Ky) G R^~^ xR, and F denotes the d-dimensional Fourier transform, 

(F(^)(K):=/ e-^^Xx)6ix, K = i^^) G R^ . 

Equation (1) has been derived in [12, 13] for three spatial dimensions. It can be 
proven in any dimension by using the inversion formula for the spherical mean 
Radon transform of [7, 9] . A related formula using the Fourier cosine transform 
instead of the Fourier transform has been obtained in [34, 15] for (i = 2, 3. 



2.2 Partial Fourier reconstruction 

The inversion formula (1) yields an exact reconstruction of /, provided that 
(Q/)(x,t) is given for all G R^. In practical applications only a partial 

(or limited view) data set is available [35, 36, 37]. In this paper we assume 
that data (Q/)(x,t) are given only for G (0,X)'^, see Figure 1, which are 
modeled by 

g{x,t) :=w,^t{x,t){Qf){x,t), (2) 

where Wcut is a smooth nonnegative cutoff function that vanishes outside (0, X)^. 
Using data (2), the function / cannot be exactly reconstructed in a stable way 
(see [38, 37]). It is therefore common to apply the exact inverse of Q to the 
partial data g and to consider the result as an approximation of the object to 
be reconstructed. More precisely, the function /''' defined by 



2Ky{Fg) lK,,sign{Ky)J\K,\^+Kl] 

{Fp){K,,Ky) := ^ , ^ ^, (3) 

signiKy)^\K,\^+K^ 
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is considered an approximation of /. The function is called partial Fourier 
reconstruction. 

Fourier reconstruction algorithms in PAI name numerical implementations 
of (3). In this paper we apply the one-dimensional nonuniform FFT to derive a 
fast and accurate algorithm for implementing (3). 



3 The Nonuniform Fast Fourier Transform 

The discrete Fourier transform of a vector g = {gn)n^Q G with respect to 
the nodes u) = (^/c)^?_7v/2 (^i^^ ^ even) is defined by 

N-1 

T[g]K):= ^e-^-^^^-/^, fe = -7V/2,...,7V/2-l. (4) 

n=0 

Direct evaluation of the N sums in (4) requires 0{N'^) operations. Using the 
classical fast Fourier transform (FFT) this effort can be reduced to 0{N log N) 
operations. However, application of the classical FFT is restricted to the case 
of evenly spaced nodes uok = k = — iV/2, . . . , N/2 — 1. 

The one-dimensional nonuniform FFT (see [18, 19, 20, 29, 21, 22, 23]) is an 
approximate but highly accurate method for evaluating (4) at arbitrary nodes 
cok, k = -N/2,...,N/2 - 1, in 0{N log N) operations. 



3.1 Derivation of the nonuniform FFT 

To derive the nonuniform FFT we closely follow the presentation of [29], which 
is based on the following lemma: 

Lemma 1 ([29, Proposition 1]). Let c > 1 and a < 7r(2c — 1). Assume that 
: R ^ R is continuous in [—a^a], vanishing outside [—a^a], and positive in 
[— TT, tt] . Then 

^^(c^-j7c)e-^^-^/% ueR,\0\<7r. (5) 



27r^((9) 

Here ^(cj) := j^e~'^^^^{0)dO denotes the one- dimensional Fourier transform 
of^. 

Proposition 2. Let c, a, ^ , and ^ he as in Lemma 1. Then, for every g = 
{gn)n=o ^ u ^R we have 

N-1 

^-ic.n2n/Ng^ ^ e-^H^-j/c)^^^ _ j/c)Gj , (6) 

n=0 jeZ 



with 

c a e-^^'^27r/(A/"c) 
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Proof. Taking = n27r/N — tt e [— 7r,7r] in (5), gives 



27r^(n27r/7V-7r) 



and therefore 



N-1 N-1 



^ 27r^^ ^ ^'^(n27r N -tt) 

n=0 n=0 jGZ ^ ' ' 

Interchanging the order of summation in the right hand side of the above equa- 
tion shows (6), (7) and concludes the proof. □ 

In the fohowing we assume that cN is an even number. Then 

/cN-l \ 

where Qn := for n > TV, is an oversampled discrete Fourier transform with the 
oversamphng factor c. Moreover we assume that ^ is concentrated around zero 
and decays rapidly away from zero. The nonuniform FFT uses the formulas (6), 
(8) to evaluate T[g] at the nodes uok- The basic steps of the algorithm are as 
follows: 

(i) Append (c — 1)A^ zeros to the vector g = {gn)n^o evaluate Gj, j = 
-Nc/2, . . . , Nc/2 - 1, in (8) with the FFT algorithm. 

(ii) Evaluate the sums in (6) approximately by using only the terms with 

— j/c\ ^ where the interpolation length K is di small positive 
parameter. 

Since ^ is assumed to decay rapidly, the truncation error in Step (ii) is small. 

The nonuniform Fourier transform is summarized in Algorithm 1 . All evalu- 
ations of ^ and 4^ are precomputed and stored. Moreover, the classical FFT is 
applied to a vector of length cN. Therefore the numerical complexity of Algo- 
rithm 1 is O(cA^logA^). Typically c = 2, in which case the numerical effort of 
the nonuniform FFT is essentially twice the effort of the one-dimensional clas- 
sical FFT applied to an input vector of the same length. See [29, Section 3] for 
an exact operation count, and a comparison between actual computation times 
of the classical and the nonuniform FFT. 

3.2 Kaiser Bessel window 

In our implementation we choose for ^ the Kaiser Bessel window^ 



K,^,._ 1 [ Io{KVa^-0^), if|^|<a, 
Io{aK) \ 0, if \e\>a. 
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Algorithm 1 Nonuniform FFT with respect to the nodes uj = (^fe)^?. 7x^/2' 
using input vector g = {gn)n^o ^ oversampUng c > 1, interpolation length 
and window function ^. 



* ^ (^(27rn/A/'-7r))^ 

function nuf f t(g, u;, c, iiT, ^, 4^) 
g ^ g/* • c/(2^) 
g ^ (g,zeros(l,(c- 1)N) 
g ^ f f t(g) 

for k = -N/2, . . . , 7V/2 - 1 do 

9k ^ Yj\j-cuJk\<cK^k,jgj 

end for 
return {gk)k 
end function 



> precomputations 



> zero-padding 
[> one-dimensional FFT 

D> interpolation 



Here /q is the modified Bessel function of order zero. The one-dimensional 
Fourier transform of ^k'b^ 

^^'^(cj) = 2sinh(a/^^^^)/(/o(ai^)/^^^^) , 

if G M \ {-K, K} and 2a/{Io{aK)) otherwise. 

The Kaiser Bessel window is a good and often used candidate for ^, since 
^KB^(^) becomes extremely small for > K. For example, with the parame- 
ters K = 3, and a = Stt, we have for uo > 

I*k'b''(^)/*k^(0)I < |*k'^ W/*k'^(0)| ^ 3 * 10-" . 

Remark 3. Take c = 1 and let ^ be the characteristic function of the interval 
[— 7r,7r]. Then ^(cj) = 27r sinc(7ra;) and (6), (7) reduce to the sine series 

N-l /N-1 \ 

^-^c.n2n/Ng^ ^ ^ g-^-(--j) ^^^^^^ ' j) [ ^ne"'^"'"/"^ , 

n=0 jez \n=0 / 

i(;/iic/i is a discretized version of Shannon's sampling formula [39, 40] 
g{u) = ^ e-^-(--^') sinc(c^ - j)g{j) , G R , 

applied to the Fourier transform of a function ^ : R ^ R that vanishes outside 
[0,27r]. 

See Figure 2 for a comparison of sine and ^kb^^ '^^^^ K = 3 and a = Stt. 
One realizes that ^^b" decays much faster than sine and is therefore much better 
suited for truncated interpolation. In fact, |^^'^(cc;)/^^'^(0)| < 3 * 10""^"^ for 
(jO >3, whereas \ sinc(a;)| < 0.01 only for uo > IOO/tt. 
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-6 -4 -2 2 4 6 '^^-^10 -5 5 10 



spatial variable 9 frequency variable w 

Figure 2: Left: Kaiser-Bessel window ^^'^(6>) and characteristic function of 
the interval [— 7r,7r]. Right: Fourier transforms ^^'^(cj) and 27rsinc(cj) in dB 
(decibel). Here dB denotes the logarithmic decay 10 logio(|(/>(^)/(/>(0)|) of some 
quantity ^(cc;). 



An error estimate for the nonuniform FFT using the Kaiser Bessel window 
is given in [29]. The result is 



27r^((9) 

^ ^ \u;-j/c\<K 



J2 ^'^^{co-j/c)e- 



ijO/c 



30 

< 



For example, taking c = 2, a = 37r and K = 3, the above error is as small as 
3* 10-^ 



4 A Fourier reconstruction Algorithm based on 
the nonuniform FFT 

In this section we apply the one-dimensional nonuniform FFT to photoacoustic 
imaging. Throughout the following we restrict our attention to two dimensions, 
noting that the general case d > 2 can be treated in an analogous manner. 

Assume that / is a smooth function that vanishes outside (0,X)^, and set 
g := Went Q/, where Wcut is as in (2). Fourier reconstruction names an imple- 
mentation of (3), that uses discrete data 

9m,n di'^^sa.mp') '^^samp) 7 (9) 

with (m, n) G {0, . . . , — 1}^ and reconstructs an approximation 

fm,n — ('^^samp? "^^samp) •> (-^0) 

with (m, n) G {0, . . . , AT - Here p is defined by (3), A' is an even number 
and Agamp •= X/N. In the Appendix we show that the sampling in (10), (9) is 
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sufficiently fine, provided that Agamp < tt/I^, where ft is the essential bandwidth 
of/. 

Discretizing (3) with the trapezoidal rule gives 

i{ln-\-km)27T/N r \ 
Jm,n I 

, N-l /N-1 \ 

= ^Y1 e-^-'-^-/^ Yl e-^'='"^-/^5„.,n , (11) 

' n—0 \m=0 / 

where 

:=sign(/)V^2^/2^ (^,0 G {-7V/2, . . . , 7V/2 - 1}^ 
One notices that the inner sums in (11), 

N-l 

~9k,n ■■= e-'""'"/'^5m.n (12) 

m=0 

can be exactly evaluated with N one-dimensional FFTs, and the outer sums 

N-l 

9kM := ^~'^'''"^^^''~9Kn (13) 

n=0 

can be approximately evaluated with TV one-dimensional nonuniform FFTs. De- 
noting the resulting approximation by Qk^i — gk{^k,i) and setting 

:= ^ , (fc, e {-N/2, ...,N/2-l}\ (14) 

^k,l 

we finally find 

N/2-1 

fn,m:=^ E 6^^^-+^-)^-/^ A, (15) 
k,l^N/2 

with the inverse two-dimensional FFT algorithm. 

The nonuniform FFT based reconstruction algorithm is summarized in Al- 
gorithm 2. Its numerical complexity can easily be estimated. Evaluating (12) 
requires NO{N log N) operations {N one-dimensional FFTs), evaluating (13) 
requires A^O(A^logA^) operations {N non-uniform FFTs), and (15) is evaluated 
with the inverse two-dimensional FFT in 0{N'^ log N) operations. Therefore 
the overall complexity of Algorithm 1 is 0{N'^ log N). 

In the next section we numerically compare Algorithm 2 with standard 
Fourier algorithms presented in the literature [41, 13], which all differ in the 
way how the sums in (13) are evaluated: 

1. Direct Fourier algorithm. Equation (13) cannot be evaluated with 
the classical FFT algorithm because the nodes u;k,i are non-equispaced. 



AT-l /N-l 

EE' 

n=0 \m=0 
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Algorithm 2 Nonuniform FFT based algorithm for calculating f = {fm,n)n m=o 



using data g = {gm,n)m n=0' oversampling factor c, interpolation length and 
window function ^. 



^ ^ (^(27rn/iV-7r))^ 

function FouRecNuf f t(g, c, K", ^, 4^) 
forn = 0, . . . , TV - 1 do 

h ^ (^m,n)m 
(^/c,n)fe ^ f ft(h) 

end for 

l^(-7V/2,...,7V/2-l) 
for k = -N/2, . . , A/'/2 - 1 do 
^ sigii(l)VA:2 +12 

h ^ nuf f t(h, u;, c, i^, ^, ^) 

{fk,i)i^2kh/u; 
end for 

f ^ {fk,l)k,l 

f ^ if f t2(f) 
return f 
end function 



> precomputations 



> one-dimensional FFT 



> nonuniform FFT 



two-dimensional inverse FFT 



The most simple way to evaluate (13) is with direct summation. Because 
there are TV^ such sums, direct Fourier reconstruction requires 0{N^) 
operations. Consequently it does not lead to a fast algorithm. However, 
since (13) is evaluated exactly, it is optimally suited to evaluate the image 
quality of reconstructions with fast Fourier algorithms. 

2. Interpolation based algorithm. A fast and simple alternative to direct 
Fourier reconstruction is as follows: Choose an oversampling factor c > 1 
and exactly evaluate 

N-1 
n=0 

at the uniformly spaced nodes uo = Agampj/c, j G {0, . . . , Nc — 1}, with 
the one-dimensional FFT algorithm. In a next step, linear interpolation 
is used to find approximate values gk^i — gk{^k,i)^ see [13]. Evaluating 
gk,i with linear interpolation requires 0{N'^) operation and therefore the 
overall numerical effort of linear interpolation based Fourier reconstruction 
is 0{N^ log TV). 

Algorithms using nearest neighbor interpolation instead of the linear one 
have the same numerical complexity and have also been applied to PAI 
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(see, e.g., [42]). Higher order polynomial interpolation has been applied 
in [43] for a cubic recording geometry. 

3. Truncated sine reconstruction. If the function ^ in Algorithm 2 is 
chosen as the characteristic function of the interval [— C7r,C7r], c > 1, then 
the nonuniform fast Fourier transform reduces to the truncated sine in- 
terpolation considered in [41]. However, due to the slow decay of sinc(c(;), 
truncation will introduce a non-negligible error in the reconstructed image 
(see Remark 3). 

The Fourier algorithms are also be compared with a numerical implementa- 
tion of the back-projection formula 

where r = {x — x'Y' + denotes the distance between the detector location 
(x',0) and the reconstruction point (x^y). Equation (16) has been obtained 
in [8] by applying the method of descent to the three-dimensional universal 
back-projection formula discovered by Xu and Wang [10]. Again (16) gives an 
exact reconstruction only if it is applied to complete data (Q/)(x,t), (x,t) G 
M^. In the numerical experiments the back-projection formula is applied to the 
partial data WcutQf^ see (2), and implemented with 0{N^) operation counts 
as described in [8, Section 3.3]. 




5 Numerical Results 

In the following we numerically compare the proposed nonuniform FFT based 
algorithm with standard Fourier algorithm and the back projection algorithm 
based on (16). 

The cutoff function Wcut is constructed by convolution of 



(Pe{x,t) 



'C,exp{-l/{e-x^ -t^f), ifx^^t^ < e, 
0, otherwise , 



with the characteristic function of [0,X]^, where e is a small parameter and Ce 
is chosen in such a way that J^2 (Pe{x^t)dxdt = 1. Typically, e is chosen as a 
"small" multiple of the sampling step size Agamp = X/N. 

In all numerical experiments we take X = 1, and TV = 512. The win- 
dow width a is chosen to be slightly smaller than 7r(2c — 1), where c is the 
oversampling factor that determines the accuracy of the Fourier reconstruction 
algorithms. 
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nearest, c = 1 H linear, c = 1 




nearest, c = 2 H linear, c = 2 

i I I 



Figure 3: Reconstruction with interpolation based Fourier algorithms. 

White corresponds to function value 1, black to function value -0.4. Top Line: 
Phantom and analytic data. Middle line: Reconstruction without oversampling 
(c = 1). Bottom line: Reconstruction with oversampling (c = 2). 



5.1 Circular shaped object 

As first case example we use a circular shaped object 

/^.^^(x) = -| (^^^ - |x-xop)^^^ if|x-xo|<a, 
^^^^ a 1 0, otherwise, 

centered at xq := (a^c^o) (see top left image in Figure 3). For such a simple 
object reconstruction artifacts can be identified very clearly. Moreover, the data 
Q/circ can be evaluated analytically (see [8, Equation (B.l)]) as 

(Q/circ)(^,^) = - Re 



.(^+-^-)-^^^H ^- + (^-ao j 
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back projection 



direct 



— L- 

truncated sine, c = 2 ■ nonuniform FFT, c = 2 




Figure 4: Improved reconstructions. Top Line: Back projection (left) and 
direct reconstruction (right). Middle Line: Truncated sine (left) and nonuni- 
form FFT based Fourier algorithm (right). Here white corresponds to function 
value 1, black to function value -0.4. Bottom Line: Difference images between 
direct and truncated sine reconstruction (left), and direct and nonuniform FFT 
based reconstruction (right). Here white (resp. black) corresponds to function 
value 0.04 (resp. -0.04). 



Here s± := {{t ± a)^ + 1(^,0) — xop)"*^/^, log(-) is the principal branch of the 
complex logarithm, and Re[2;] denotes the real part of complex number z. The 
reconstruction results are depicted in Figures 3 and 4. Table 1 and Figure 5 
compare run times with the relative £^-error 

/ \ 1/2 

||f - ft 11^2 y^m,nifm,n - /m,n)^j 
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Figure 5: Reconstruction time versus error. The points on the graphs 
belong to runtimes and errors for reconstruction with oversamphng factors c G 
{1,2,4,8,16,32}. 



were f ^ = (/^ ^) denotes the discrete image obtained by direct Fourier recon- 
struction. Run times were measured for Matlab implementations on a personal 
computer with 2.4 GHz Athlon processor. 

In order to demonstrate the stability of the Fourier algorithms, we also 
performed reconstructions from noisy data, where Gaussian noise was added 
with a variance equal to 20% of the maximal data value. The reconstruction 
results are depicted in Figure 7. 





c 


^^-error 


runtime (sec) 


back projection 






88.9 


direct reconstruction 






54.1 


nearest neighbor 


1 


0.75 


0.65 


nearest neighbor 


2 


0.40 


0.85 


linear 


1 


0.65 


0.8 


linear 


2 


0.21 


0.95 


Truncates sine 


2 


0.04 


1.6 


Kaiser Bessel 


2 


0.006 


1.6 



Table 1: Run times and error of different reconstruction methods. 
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Figure 6: Noisy data used in Figure 7 



5.2 Shepp— Logan phantom 

In the next example we consider the Shepp-Logan phantom /phant, which is 
shown in top left image in Figure 8. The data were calculated numerically by 
implementing d'Alemberts formula [44], 

(Q/phant)(^,0,t) =dt _ dr 

with 

(M/phant)(^,0,r) — (p /phant((3^,0) + m) dcr 

denoting the spherical mean Radon transform of fdrc The reconstruction re- 
sults from simulated data are depicted in Figure 8. 

5.3 Discussion 

We emphasize that none of the above Fourier algorithms are designed to cal- 
culate an approximation of / but an approximation to the partial Fourier re- 
construction defined in (3). Therefore even in the direct reconstruction (top 
right image in Figure 4) and in the back projection reconstruction one can see 
some blurred boundaries in the reconstruction. Such artifacts are expected using 
limited view data (2), see [38, 37]. 

The results of interpolation based reconstruction without oversampling (c = 
1) are quite useless. The reconstructions are significantly improved by using a 
larger oversampling factor c. However, even then, the results never reach the 
quality of the nonuniform FFT based reconstruction. Moreover, the numerical 
effort of linear interpolation based reconstruction is proportional to the over- 
sampling factor, which prohibits the use of "very large" values for c (see Figure 
5). In the reconstruction with c = 2 (bottom line in Figure 3 and middle line 
in Figure 8) artifacts are still clearly visible. 

The images in the middle line of Figure 4 suggest that truncated sine and 
nonuniform FFT based reconstruction seem to perform quite similar. However, 
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back projection 



• 



100 200 300 400 500 



nonuniform FFT, c = 2 



ill I exact I" 

A ^ 

||#*'^ 



100 200 300 400 500 

Figure 7: Reconstruction from noisy data. Top line: Nearest neighbor 
interpolation based reconstruction (c = 2). Second line: Linear interpolation 
based reconstruction (c = 2). Third line: Reconstruction with back projection 
algorithm. Bottom line: Reconstruction with nonuniform FFT algorithm (c = 
2). The horizontal profiles on the right are taken at {x = xq}. 
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Figure 8: Reconstruction of Shepp Logan phantom. Top line: Phantom 
and simulated data. Second line : Interpolation based reconstruction. Bot- 
tom line: Reconstruction with back projection algorithm (left) and proposed 
nonuniform FFT algorithm (right). 

the differences to the direct Fourier reconstructions, shown in the bottom line 
in Figure 4, demonstrate the higher accuracy of the nonuniform FFT based 
algorithm. 

The results in Figure 7 show that all applied reconstruction algorithms are 
quite stable with respect to data perturbation. In particular, the filtered back 
projection algorithm produces images with the highest signal to noise ratio. 
However, only at the cost of a nearly 100 times longer computation time (see 
Table 1). 
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6 Conclusion 



We presented a novel fast Fourier reconstruction algorithm for photoacoustic 
imaging using a limited planar detector array. The proposed algorithm is based 
on the nonuniform FFT. Theoretical investigation as well as numerical simu- 
lations show that our algorithm produces better images than existing Fourier 
algorithms with a similar numerical complexity. Moreover the proposed algo- 
rithm has been shown to be stable against data perturbations. 
[Sampling and Resolution] 

Let / be smooth function that vanishes outside (0,X)^, and define 
by (2), (3). We further assume that F Wcut is concentrated around zero and, 
that / is essentially bandlimited with essential bandwidth in the sense that 
(F/)(K) is negligible for |K| > Note that since / has bounded support, F f 
cannot vanish exactly on {|K| > Q}. 

• Sampling of g. Equation (1) implies that 

(F g){K,,u;) = (F w,^t) * (F Q /) {K,,u;) , (17) 

with 

(FQ/)(i^.,c.) = .\ . . , ,^ ^ 

sign(cj)Y/cj^ - \Kx\^ 

if |cc;| > \Kx\'^ and {F Q f){Kx,uj) = otherwise. The assumption that / 
has essential bandwidth Q and equation (17) imply that (F g){Kx^uj) is 
negligible outside the set 

{{Kx,co):\Kx\<\co\<n}c{-n,ny. 

Now Shannon's sampling theorem [39, 40] states that g is sufficiently 
fine sampled if the step size in x and in t satisfies the Nyquist condition 

• Sampling of /'''. Similar considerations as above again show that /''' 
is essentially bandlimited with essential bandwidth Q. Shannon's sam- 
pling theorem implies that can be reliable reconstructed from discrete 
samples taken with step size Agamp < tt/Q. 

If / has essential bandwidth larger than the function g has to be filtered 
with a low pass-filter before sampling. Otherwise, sampling introduces aliasing 
artifacts in the reconstructed image. 

Theoretically, the resolution (at least of the visible parts) can be increased ad 
infinity by simply decreasing the sampling size Agamp- In practical applications, 
several other factors such as the bandwidth of the ultrasound detection system 
limit the bandwidth of the data, and therefore the resolution of reconstructed 
images [33]. This, however, also guarantees that in practice a moderate sampling 
step size Agamp gives correct sampling without aliasing. 
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